# Install 'pacman' once if you do not have it (uncomment next line and run only once)
# install.packages("pacman")
# Use pacman::p_load() to install (if needed) and then load the listed packages automatically
pacman::p_load(
sf, # modern vector geospatial data handling (Simple Features)
terra, # modern raster + vector geospatial data handling
spatstat, # meta‑package: tools for spatial point pattern analysis (SPPA)
tmap, # cartographic and interactive mapping
rvest, # web‑scraping helper used here to parse HTML inside KML attributes
tidyverse # data wrangling helpers (dplyr, purrr, stringr, readr, etc.)
)Hands-on_Ex02
4. 1st Order Spatial Point Patterns Analysis Methods
4.1 Overview
The specific questions we would like to answer are as follows:
- Are the childcare centres in Singapore randomly distributed throughout the country?
- If the answer is not, then the next logical question is where are the locations with higher concentration of childcare centres?
4.3 Installing and Loading the R packagesn
4.4 Importing and Wrangling Geospatial Data Sets
Import the Master Plan 2019 Subzone (No Sea) polygons as sf and set projection to EPSG:3414 (SVY21 / Singapore TM)
# Read the subzone boundary KML into an sf object named 'mpsz_sf'
mpsz_sf <- sf::st_read("data/geospatial/MasterPlan2019SubzoneBoundaryNoSeaKML.kml") %>%
sf::st_zm(drop = TRUE, what = "ZM") %>% # remove Z (elevation) and M (measure) dimensions to keep 2D only
sf::st_transform(crs = 3414) Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source
`/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/Hands-on_Ex02/Data/geospatial/MasterPlan2019SubzoneBoundaryNoSeaKML.kml'
using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY, XYZ
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84
Why remove Z/M? Our analysis is purely 2D on a planar map. Extra dimensions are unnecessary and can break some operations.
Why EPSG:3414? This is the standard projected system for Singapore. Distances/areas become meaningful (meters) instead of degrees.
Extract REGION_N, PLN_AREA_N, SUBZONE_N, SUBZONE_C from the Description field (HTML inside KML) (as in the guidance helper function)
# Define a small helper function to pull one named item from the HTML table stored in the 'Description' field
extract_kml_field <- function(html_text, field_name) {
if (is.na(html_text) || html_text == "") return(NA_character_) # if Description is empty, return NA immediately
page <- rvest::read_html(html_text) # parse the HTML string
rows <- rvest::html_elements(page, "tr") # get all table rows (<tr>)
# Find the row whose header cell (<th>) equals the field_name, then take the value in its <td>
value <- rows %>%
purrr::keep(~ rvest::html_text2(rvest::html_element(.x, "th")) == field_name) %>%
rvest::html_element("td") %>%
rvest::html_text2()
if (length(value) == 0) NA_character_ else value # if not found, return NA
}# Apply the function to build clean attributes from 'Description'
mpsz_sf <- mpsz_sf %>%
dplyr::mutate(
REGION_N = purrr::map_chr(Description, extract_kml_field, "REGION_N"), # region name
PLN_AREA_N = purrr::map_chr(Description, extract_kml_field, "PLN_AREA_N"), # planning area name
SUBZONE_N = purrr::map_chr(Description, extract_kml_field, "SUBZONE_N"), # subzone name
SUBZONE_C = purrr::map_chr(Description, extract_kml_field, "SUBZONE_C") # subzone code
) %>%
dplyr::select(-Name, -Description) %>% # drop original noisy fields
dplyr::relocate(geometry, .after = dplyr::last_col()) # move geometry column to the end for readabilityCheckpoint: You now have clean polygon attributes and the proper CRS.
(Optional) Filter out offshore/irrelevant areas exactly like in the slides
# Create a cleaned version that excludes southern group, western islands, and north‑eastern islands (as per slides)
mpsz_cl <- mpsz_sf %>%
dplyr::filter(
SUBZONE_N != "SOUTHERN GROUP",
PLN_AREA_N != "WESTERN ISLANDS",
PLN_AREA_N != "NORTH-EASTERN ISLANDS"
)
# Save a small RDS for repeatability (optional)
# readr::write_rds(mpsz_cl, "chap04/data/mpsz_cl.rds")Import Child Care Services points as sf, drop Z/M, and project to EPSG:3414
# Read childcare locations KML into an sf object named 'childcare_sf'
childcare_sf <- sf::st_read("data/geospatial/ChildCareServices.kml") %>%
sf::st_zm(drop = TRUE, what = "ZM") %>% # keep 2D only
sf::st_transform(crs = 3414) Reading layer `CHILDCARE' from data source
`/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/Hands-on_Ex02/Data/geospatial/ChildCareServices.kml'
using driver `KML'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84
4.4.1 Mapping the Geospatial Data Sets
# Quick mapping of both layers to verify alignment (childcare points inside planning subzones)
tmap::tmap_mode("plot") # ensure static plotting modeℹ tmap mode set to "plot".
# Plot subzones with childcare points overlaid
tmap::tm_shape(mpsz_cl) +
tmap::tm_polygons(col = "grey85", border.col = "black") +
tmap::tm_shape(childcare_sf) +
tmap::tm_dots(size = 0.2, col = "red")
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
(instead of 'col'), and 'col' for the outlines (instead of 'border.col').

What to check: All childcare centres (black dots) fall neatly inside Singapore’s subzones (grey polygons). If not, check CRS and transformations.
Interactive mapping with tmap + leaflet
# Switch tmap into interactive mode (internally uses the leaflet library)
tmap::tmap_mode("view") # enable pan/zoom and feature popupsℹ tmap mode set to "view".
# Draw an interactive point map on top of a web basemap
# - Click a point to see its attributes
# - Use the layer control to change basemaps (e.g., ESRI.WorldGrayCanvas, OpenStreetMap)
tmap::tm_shape(childcare_sf) + # select the childcare sf layer
tmap::tm_dots() # render as clickable pointsRegistered S3 method overwritten by 'jsonify':
method from
print.json jsonlite
# IMPORTANT: After exploring, switch back to static plotting for the rest of the workflow
# This avoids keeping too many active web-map connections during knitting/deployment
tmap::tmap_mode("plot") # restore static plotting modeℹ tmap mode set to "plot".
# Switch to interactive mode for web-map style exploration
# tmap::tmap_mode("view")
#
# # Create an interactive map where you can zoom and click points
# tmap::tm_shape(mpsz_cl) +
# tmap::tm_polygons(col = "lightgrey", border.col = "white") +
# tmap::tm_shape(childcare_sf) +
# tmap::tm_dots(col = "blue", size = 0.5)
#
# # Always switch back to static plotting when done
# tmap::tmap_mode("plot")Notes: In interactive mode you can zoom, pan, and click childcare points to view their attributes. Remember to switch back to plot mode before continuing with static maps.
4.5 Geospatial Data Wrangling
4.5.1 Convert sf Point Data to ppp (planar point pattern)
# Convert the sf point layer to spatstat's planar point pattern (ppp) object
childcare_ppp <- spatstat.geom::as.ppp(childcare_sf) # now suitable for SPPA functions# Verify the class
class(childcare_ppp) # expect: "ppp"[1] "ppp"
# Peek at summary to understand the point pattern and its window
summary(childcare_ppp) # quick glance (number of points, intensity, window size)Marked planar point pattern: 1925 points
Average intensity 2.417323e-06 points per square unit
Coordinates are given to 11 decimal places
Mark variables: Name, Description
Summary:
Name Description
Length:1925 Length:1925
Class :character Class :character
Mode :character Mode :character
Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
(33590 x 23700 units)
Window area = 796335000 square units
4.5.2 Create an Owin (observation window) from the Singapore Boundary
# Convert the cleaned subzones into a single observation window (owin) for clipping/analysis
sg_owin <- spatstat.geom::as.owin(mpsz_cl)# Verify the class
class(sg_owin) # expect: "owin"[1] "owin"
# Visualize the boundary to confirm it looks right
plot(sg_owin, main = "sg_owin — Singapore Observation Window")
4.3 Combine the ppp with the owin (keep only points inside Singapore)
# Clip the point pattern to the Singapore window to exclude any points outside the boundary
childcareSG_ppp <- childcare_ppp[sg_owin]# Confirm the combined object (ppp with polygon window information)
childcareSG_ppp # should report a polygonal window and the point countMarked planar point pattern: 1925 points
Mark variables: Name, Description
window: polygonal boundary
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units
4.6 Clark-Evan Test for Nearest Neighbour Analysis
Hypotheses:
\(H_0\): Childcare centres are randomly distributed (Complete Spatial Randomness, CSR).
\(H_1\): Childcare centres are not randomly distributed (clustered or regular).
The 95% confident interval will be used.
The Clark–Evans test returns an index R:
- \(R < 1 → clustered\).
- \(R = 1 → random\)
- \(R > 1 → regular\)
4.6.1 Clark–Evans without CSR Simulation (Z-test)
# Run Clark–Evans using observed pattern only (Z-test)
ce_noCSR <- spatstat.explore::clarkevans.test(
X = childcareSG_ppp, # ppp object of childcare centres
correction = "none", # no edge correction (per slides)
alternative = "clustered" # one-sided: test for clustering (R < 1)
)
# Show results
ce_noCSR
Clark-Evans test
No edge correction
Z-test
data: childcareSG_ppp
R = 0.53532, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
- R = 0.53532 (< 1) → indicates clustering.
- p-value < 2.2e-16 (< 0.05): → reject H₀.
- Statistical conclusion: Distribution is not random; strong clustering exists.
- Business communication: Childcare centres are concentrated in specific neighbourhoods. This reflects demand-driven placement but may leave some areas underserved. Authorities should consider equity in future allocations.
4.6.2 Clark–Evans with CSR via Monte‑Carlo (nsim = 99)
# Monte‑Carlo Clark–Evans: compare observed R to R from CSR simulations (nsim = 99)
ce_MC <- spatstat.explore::clarkevans.test( # compute simulated p‑value under CSR
X = childcareSG_ppp, # same ppp as above (metres)
correction = "none", # keep consistent with slides
alternative = "clustered", # one‑sided for clustering (R < 1)
method = "MonteCarlo", # enable Monte‑Carlo simulation mode
nsim = 99 # number of CSR replicates as per slides
) # end clarkevans.test call
ce_MC
Clark-Evans test
No edge correction
Monte Carlo test based on 99 simulations of CSR with fixed n
data: childcareSG_ppp
R = 0.53532, p-value = 0.01
alternative hypothesis: clustered (R < 1)
- R = 0.53532 (< 1): → clustering again.
- p-value = 0.01 (< 0.05) → reject H₀ at 95% CI.
- Statistical conclusion: Distribution is not random; clustering is statistically significant even under CSR simulations.
- Business communication: Confirms earlier finding with stronger robustness. Clustering is not due to chance — it is systematic. Policymakers should expand childcare access in low-density areas to reduce inequality in service availability.
4.7 Kernel Density Estimation (KDE) Method
Goal: Turn points into a smooth surface to reveal hotspots. Bandwidth controls smoothness; kernel controls the spread shape.
4.7.1 Working with automatic bandwidth selection method
kde_SG_diggle <- density( # Create KDE surface and save to object 'kde_SG_diggle'
childcareSG_ppp, # Input: childcare centres dataset in 'ppp' format
sigma = bw.diggle, # Bandwidth: automatic smoothing radius (Diggle’s selector)
edge = TRUE, # Apply edge correction to fix boundary underestimation
kernel = "gaussian" # Kernel type: Gaussian (bell-shaped smoothing function)
) # End of density() call# Plot the kernel density estimation surface for childcare centres
plot(kde_SG_diggle) 
summary(kde_SG_diggle)real-valued pixel image
128 x 128 pixel array (ny, nx)
enclosing rectangle: [2667.538, 55941.94] x [21448.47, 50256.33] units
dimensions of each pixel: 416 x 225.0614 units
Image is defined on a subset of the rectangular grid
Subset area = 669941961.12249 square units
Subset area fraction = 0.437
Pixel values (inside window):
range = [-5.824417e-21, 3.063698e-05]
integral = 1927.788
mean = 2.877545e-06
# Calculate optimal bandwidth (smoothing parameter) using Diggle’s method
bw <- bw.diggle(childcareSG_ppp)
bw # Display the selected bandwidth value sigma
295.9712
4.7.2 Rescalling KDE values
childcareSG_ppp_km <- rescale.ppp( # Rescale the point pattern dataset from metres to kilometres
childcareSG_ppp, 1000, "km") # Display the selected bandwidth value kde_childcareSG_km <- density( # Compute kernel density estimation on rescaled dataset
childcareSG_ppp_km, # Input point pattern in kilometres
sigma = bw.diggle, # Use Diggle’s optimal bandwidth (sigma)
edge = TRUE, # Apply edge correction to fix boundary bias
kernel = "gaussian" # Use Gaussian kernel for smoothing
)# Plot KDE surface for childcare centres (units in km)
plot(kde_childcareSG_km)
4.7.3 Working with different automatic badwidth methods
# Compute bandwidth using Cronie & van Lieshout method
bw.CvL(childcareSG_ppp_km) sigma
4.357209
# Compute bandwidth using Scott’s rule-of-thumb method
bw.scott(childcareSG_ppp_km) sigma.x sigma.y
2.159749 1.396455
# Compute bandwidth using likelihood cross-validation (PPL)
bw.ppl(childcareSG_ppp_km) sigma
0.378997
# Re-compute bandwidth using Diggle’s method (on km scale)
bw.diggle(childcareSG_ppp_km) sigma
0.2959712
kde_childcareSG_ppl <- density( # KDE using PPL bandwidth for smoothing
childcareSG_ppp_km, # Input dataset in kilometres
sigma = bw.ppl, # Use bandwidth chosen by PPL
edge = TRUE, # Apply edge correction
kernel = "gaussian" # Gaussian smoothing kernel
)
par(mfrow=c(1,2)) # Set plot area into 1 row, 2 columns for comparison
plot(kde_childcareSG_km, main = "bw.diggle") # Plot KDE using Diggle bandwidth
plot(kde_childcareSG_ppl, main = "bw.ppl") # Plot KDE using PPL bandwidth
4.7.4 Working with different kernel methods
par(mfrow=c(2,2)) # Divide plotting window into 2 rows × 2 columns
plot(density(childcareSG_ppp_km, # KDE with Gaussian kernel
sigma=0.2959712, edge=TRUE,
kernel="gaussian"),
main="Gaussian") # Title for Gaussian plot
plot(density(childcareSG_ppp_km, # KDE with Epanechnikov kernel
sigma=0.2959712, edge=TRUE,
kernel="epanechnikov"),
main="Epanechnikov") # Title for Epanechnikov plot
plot(density(childcareSG_ppp_km, # KDE with Quartic kernel
sigma=0.2959712, edge=TRUE,
kernel="quartic"),
main="Quartic") # Title for Quartic plot
plot(density(childcareSG_ppp_km, # KDE with Disc kernel
sigma=0.2959712, edge=TRUE,
kernel="disc"),
main="Disc") # Title for Disc plot
4.8 Fixed and Adaptive KDE
4.8.1 Computing KDE by using fixed bandwidth
kde_childcareSG_fb <- density( # KDE using fixed bandwidth method
childcareSG_ppp_km, # Input childcare centres dataset (in km)
sigma = 0.6, # Fixed bandwidth = 0.6 km (600 metres)
edge = TRUE, # Apply edge correction at boundaries
kernel = "gaussian" # Gaussian smoothing kernel
)
plot(kde_childcareSG_fb) # Plot the fixed bandwidth KDE result
4.8.2 Computing KDE by using adaptive bandwidth
kde_childcareSG_ab <- adaptive.density( # KDE using adaptive bandwidth method
childcareSG_ppp_km, # Input childcare centres dataset (in km)
method="kernel" # Kernel smoothing method
)
plot(kde_childcareSG_ab) # Plot the adaptive bandwidth KDE result
par(mfrow=c(1,2)) # Set plotting window to 1 row × 2 columns
plot(kde_childcareSG_fb, main = "Fixed bandwidth") # Show fixed bandwidth KDE
plot(kde_childcareSG_ab, main = "Adaptive bandwidth") # Show adaptive bandwidth KDE
4.9 Plotting cartographic quality KDE map
4.9.1 Converting gridded output into raster
# Convert KDE (im class) into SpatRaster object
kde_childcareSG_bw_terra <- rast(kde_childcareSG_km)# Check the class of the raster object (should be "SpatRaster")
class(kde_childcareSG_bw_terra)[1] "SpatRaster"
attr(,"package")
[1] "terra"
# Print raster properties: resolution, extent, units, etc.
kde_childcareSG_bw_terraclass : SpatRaster
size : 128, 128, 1 (nrow, ncol, nlyr)
resolution : 0.4162063, 0.2250614 (x, y)
extent : 2.667538, 55.94194, 21.44847, 50.25633 (xmin, xmax, ymin, ymax)
coord. ref. :
source(s) : memory
name : lyr.1
min value : -4.314789e-15
max value : 3.063698e+01
unit : km
4.9.2 Assigning projection systems
# Assign projection system SVY21 / Singapore TM (EPSG:3414)
crs(kde_childcareSG_bw_terra) <- "EPSG:3414"# Re-check raster details, now with CRS applied
kde_childcareSG_bw_terraclass : SpatRaster
size : 128, 128, 1 (nrow, ncol, nlyr)
resolution : 0.4162063, 0.2250614 (x, y)
extent : 2.667538, 55.94194, 21.44847, 50.25633 (xmin, xmax, ymin, ymax)
coord. ref. : SVY21 / Singapore TM (EPSG:3414)
source(s) : memory
name : lyr.1
min value : -4.314789e-15
max value : 3.063698e+01
unit : km
4.9.3 Plotting KDE map with tmap
tm_shape(kde_childcareSG_bw_terra) + # Define raster object to be plotted
tm_raster(col.scale = # Set raster colour scheme
tm_scale_continuous(values="viridis"),
col.legend = tm_legend( # Add legend for density values
title = "Density values", # Legend title
title.size = 0.7, # Legend title text size
text.size = 0.7), # Legend labels text size
bg.color = "white", # Background colour of map
bg.alpha = 0.7, # Transparency of background
position = tm_pos_in("right","bottom"), # Place legend bottom-right
frame = TRUE) + # Draw frame around raster
tm_graticules(labels.size = 0.7) + # Add graticule grid with label size 0.7
tm_compass() + # Add compass to map
tm_layout(scale = 1.0) # Set layout scale[tm_raster()] Arguments `bg.color`, `bg.alpha`, `position`, and `frame`
unknown.
[plot mode] fit legend/component: Some legend items or map compoments do not
fit well, and are therefore rescaled.
ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

4.10 First Order SPPA at the Planning Subzone Level
4.10.1 Geospatial data wrangling
4.10.1.1 Extracting study area
pg <- mpsz_cl %>% # Create dataset 'pg' by filtering master plan polygons
filter(PLN_AREA_N == "PUNGGOL") # Keep only polygons where planning area = Punggol
tm <- mpsz_cl %>% # Create dataset 'tm'
filter(PLN_AREA_N == "TAMPINES") # Keep only polygons where planning area = Tampines
ck <- mpsz_cl %>% # Create dataset 'ck'
filter(PLN_AREA_N == "CHOA CHU KANG") # Keep only polygons where planning area = Choa Chu Kang
jw <- mpsz_cl %>% # Create dataset 'jw'
filter(PLN_AREA_N == "JURONG WEST") # Keep only polygons where planning area = Jurong Westpar(mfrow=c(2,2)) # Arrange plotting area into 2 rows × 2 columns
plot(st_geometry(pg), main="Punggol") # Plot polygon boundary of Punggol planning area
plot(st_geometry(tm), main="Tampines") # Plot polygon boundary of Tampines planning area
plot(st_geometry(ck), main="Choa Chu Kang") # Plot polygon boundary of Choa Chu Kang planning area
plot(st_geometry(jw), main="Jurong West") # Plot polygon boundary of Jurong West planning area
4.10.1.2 Creating owin object
pg_owin <- as.owin(pg) # Convert Punggol polygon(s) to an 'owin' window (study region)
tm_owin <- as.owin(tm) # Convert Tampines polygon(s) to 'owin'
ck_owin <- as.owin(ck) # Convert Choa Chu Kang polygon(s) to 'owin'
jw_owin <- as.owin(jw) # Convert Jurong West polygon(s) to 'owin'4.10.1.3 Combining point events object and owin object
childcare_pg_ppp <- childcare_ppp[pg_owin] # Clip national points to Punggol window (points inside only)
childcare_tm_ppp <- childcare_ppp[tm_owin] # Clip points to Tampines window
childcare_ck_ppp <- childcare_ppp[ck_owin] # Clip points to Choa Chu Kang window
childcare_jw_ppp <- childcare_ppp[jw_owin] # Clip points to Jurong West windowchildcare_pg_ppp_km <- rescale.ppp( # Rescale Punggol points from metres → kilometres
childcare_pg_ppp, 1000, "km" # divide coords by 1000; label new unit as "km"
)
childcare_tm_ppp_km <- rescale.ppp( # Rescale Tampines points to kilometres
childcare_tm_ppp, 1000, "km"
)
childcare_ck_ppp_km <- rescale.ppp( # Rescale CCK points to kilometres
childcare_ck_ppp, 1000, "km"
)
childcare_jw_ppp_km <- rescale.ppp( # Rescale Jurong West points to kilometres
childcare_jw_ppp, 1000, "km"
)par(mfrow = c(2,2)) # Arrange plotting area into a 2×2 grid
plot(unmark(childcare_pg_ppp_km), # Plot Punggol points (unmark = hide text marks)
main = "Punggol") # Panel title
plot(unmark(childcare_tm_ppp_km), # Plot Tampines points
main = "Tampines")
plot(unmark(childcare_ck_ppp_km), # Plot Choa Chu Kang points
main = "Choa Chu Kang")
plot(unmark(childcare_jw_ppp_km), # Plot Jurong West points
main = "Jurong West")
par(mfrow = c(1,1)) # Reset plotting layout back to single panel4.10.2 Clark and Evans Test
clarkevans.test(childcare_ck_ppp, # Clark–Evans test for the Choa Chu Kang point pattern
correction = "none", # No edge correction (consistent with slides)
clipregion = NULL, # Use the pattern's own observation window
alternative = c("two.sided"), # Two-sided hypothesis (clustered or regular)
nsim = 999 # 999 CSR simulations for p-value
) # End of clarkevans.test call
Clark-Evans test
No edge correction
Z-test
data: childcare_ck_ppp
R = 0.84097, p-value = 0.008866
alternative hypothesis: two-sided
4.10.2.2 Tampines planning area
clarkevans.test(childcare_tm_ppp, # Clark–Evans test for the Tampines point pattern
correction = "none", # No edge correction (to match slides)
clipregion = NULL, # Do not clip by an additional region (use pattern's window)
alternative = c("two.sided"), # Two-sided test: allow clustering (R<1) or regularity (R>1)
nsim = 999 # Monte-Carlo CSR simulations (999 replicates)
) # End of clarkevans.test call
Clark-Evans test
No edge correction
Z-test
data: childcare_tm_ppp
R = 0.66817, p-value = 6.58e-12
alternative hypothesis: two-sided
4.10.3 Computing KDE surfaces by planning area
par(mfrow = c(2,2)) # Arrange plotting window into 2 rows × 2 columns
plot(density(childcare_pg_ppp_km, # KDE for Punggol (data already rescaled to km)
sigma = bw.diggle, # Use Diggle’s automatic bandwidth selector
edge = TRUE, # Apply edge correction (reduces boundary bias)
kernel = "gaussian"), # Gaussian kernel (smooth bell-shaped influence)
main = "Punggol") # Panel title
plot(density(childcare_tm_ppp_km, # KDE for Tampines
sigma = bw.diggle, # Same bandwidth rule for comparability
edge = TRUE, # Edge correction on
kernel = "gaussian"), # Gaussian kernel
main = "Tampines") # Panel title
plot(density(childcare_ck_ppp_km, # KDE for Choa Chu Kang
sigma = bw.diggle, # Diggle bandwidth (km)
edge = TRUE, # Edge correction on
kernel = "gaussian"), # Gaussian kernel
main = "Choa Chu Kang") # Panel title
plot(density(childcare_jw_ppp_km, # KDE for Jurong West
sigma = bw.diggle, # Diggle bandwidth (km)
edge = TRUE, # Edge correction on
kernel = "gaussian"), # Gaussian kernel
main = "Jurong West") # Panel title
5. 2nd Order Spatial Point Patterns Analysis Methods
5.5 Second-order Spatial Point Patterns Analysis
First-order asks “where are points denser?”; second-order asks “how do points interact with each other across distance?”
We’ll use four classical functions: - G(r) — nearest-neighbour distribution (from each event to its nearest event). - F(r) — empty-space distribution (from random locations to nearest event). - K(r) — accumulates neighbours within radius r; higher than CSR suggests clustering. - L(r) = √(K(r)/π) — variance-stabilised K; plot L(r) − r to read deviations easily.
For each function your Prof shows estimation and a Monte Carlo CSR test with envelopes.
5.6 Analysing Spatial Point Process Using G-Function
5.6.1 Choa Chu Kang planning area
5.6.1.1 Computing G-function estimation
set.seed(1234) # fix the random seed so plots are reproducible for students
G_CK <- Gest(childcare_ck_ppp, correction = "border") # estimate nearest-neighbour CDF G(r) with border edge correction
plot(G_CK, xlim = c(0, 500)) # plot G(r) up to 500 m; dashed line shows CSR (Poisson) reference
5.6.1.2 Performing Complete Spatial Randomness (CSR) Test — Monte Carlo
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
\(H_0\) = The distribution of childcare services at Choa Chu Kang are randomly distributed.
\(H_1\) = The distribution of childcare services at Choa Chu Kang are not randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.
Monte Carlo test with G-fucntion
# run 999 CSR simulations; compute simulation envelopes for G(r)
G_CK.csr <- envelope(childcare_ck_ppp, Gest, nsim = 999) Generating 999 simulations of CSR ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.
Done.
# plot observed G(r), CSR expectation, and 999-sim envelopes
plot(G_CK.csr) 
If the black observed curve is mostly above the grey envelope → points are clustered at those distances. If it is below → regular/repulsive (inhibition). Inside the envelope → not significantly different from CSR at that scale.
5.6.2 Tampines planning area
5.6.2.1 Computing G-function estimation
G_tm <- Gest(childcare_tm_ppp, correction = "best") # estimate G(r) using 'best' automatic edge correction
plot(G_tm) # plot observed G(r) vs theoretical CSR G(r)
5.6.2.2 Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
\(H_o\) = The distribution of childcare services at Tampines are randomly distributed.
\(H_1\) = The distribution of childcare services at Tampines are not randomly distributed.
The null hypothesis will be rejected is p-value is smaller than alpha value of 0.001.
The code chunk below is used to perform the hypothesis testing.
G_tm.csr <- envelope(childcare_tm_ppp, Gest, # simulate CSR envelopes for G(r) at Tampines
correction = "all", # request all supported edge corrections inside envelope calc
nsim = 999) # use 999 simulations as in Prof’s slideGenerating 999 simulations of CSR ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.
Done.
plot(G_tm.csr) # display observed curve, CSR, and envelopes
Decision rule (as per slides): Reject \(H_0\): CSR if the observed curve exits the envelope and the associated p-value < 0.001 (α = 0.001).
5.7 Analysing Spatial Point Process Using F-Function
5.7.1 Choa Chu Kang planning area
5.7.1.1 Computing F-function estimation
F_CK <- Fest(childcare_ck_ppp) # estimate empty-space CDF F(r): distance from random locations to nearest facility
plot(F_CK) # plot multiple estimators and the CSR reference curve F_pois(r)
5.7.2 Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
\(H_0\) = The distribution of childcare services at Choa Chu Kang are randomly distributed.
\(H_1\) = The distribution of childcare services at Choa Chu Kang are not randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.
Monte Carlo test with F-fucntion
F_CK.csr <- envelope(childcare_ck_ppp, Fest, nsim = 999) # simulate 999 CSR patterns; compute F(r) envelopesGenerating 999 simulations of CSR ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.
Done.
plot(F_CK.csr)
- F(r) above CSR → space tends to be closer to facilities than CSR (suggests clustering).
- F(r) below CSR → larger gaps than expected (suggests inhibition).
5.7.3 Tampines planning area
****5.7.3.1 Computing F-function estimation****
F_tm = Fest(childcare_tm_ppp, correction = "best")
plot(F_tm)
5.7.3.2 Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
Ho = The distribution of childcare services at Tampines are randomly distributed.
H1= The distribution of childcare services at Tampines are not randomly distributed.
The null hypothesis will be rejected is p-value is smaller than alpha value of 0.001.
The code chunk below is used to perform the hypothesis testing.
F_tm.csr <- envelope(childcare_tm_ppp, Fest, correction = "all", nsim = 999)Generating 999 simulations of CSR ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.
Done.
plot(F_tm.csr)
5.8 Analysing Spatial Point Process Using K-Function
5.8.1 Choa Chu Kang planning area
5.8.1.1 Computing K-fucntion estimate
# estimate K(r) using Ripley (isotropic) edge correction
K_ck <- Kest(childcare_ck_ppp, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(m)")
5.8.1.2 Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.
H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.
K_tm.csr <- envelope(childcare_tm_ppp, Kest, # build CSR envelopes for K(r) in Tampines
nsim = 99, # 99 simulations exactly as shown
rank = 1, # use rank-1 global envelope
glocal = TRUE) # enable global+local envelope calculationGenerating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
99.
Done.
plot(K_tm.csr, . - r ~ r, # plot (K_hat - r) vs r using the formula used by Prof
xlab = "d", ylab = "K(d)-r", # match labels from the screenshot
xlim = c(0, 500)) # limit the x-axis to 0–500 m as in the slide
- Under CSR, K(r) = πr² and L(r) − r = 0.
- Curve above CSR → clustering; below → inhibition; inside envelopes → not significantly different from CSR.
5.8.2 Tampines planning area
5.8.2.1 Computing K-fucntion estimation
K_tm = Kest(childcare_tm_ppp, correction = "Ripley")
plot(K_tm, . -r ~ r,
ylab= "K(d)-r", xlab = "d(m)",
xlim=c(0,1000))
5.8.2.2 Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
\(H_0\) = The distribution of childcare services at Tampines are randomly distributed.
\(H_1\) = The distribution of childcare services at Tampines are not randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.
The code chunk below is used to perform the hypothesis testing.
K_tm.csr <- envelope(childcare_tm_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
99.
Done.
plot(K_tm.csr, . - r ~ r,
xlab="d", ylab="K(d)-r", xlim=c(0,500))
5.9 Analysing Spatial Point Process Using L-Function
5.9.1 Choa Chu Kang planning area
5.9.1.1 Computing L Fucntion estimation
L_ck <- Lest(childcare_ck_ppp, correction = "Ripley") # estimate L(r) using Ripley edge correction
plot(L_ck, . - r ~ r, # plot (L_hat - r) vs r to centre CSR at zero
ylab = "L(d)-r", xlab = "d(m)") # match axis labels exactly as shown
5.9.1.2 Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
\(H_0\) = The distribution of childcare services at Choa Chu Kang are randomly distributed.
\(H_1\) = The distribution of childcare services at Choa Chu Kang are not randomly distributed.
The null hypothesis will be rejected if p-value if smaller than alpha value of 0.001.
The code chunk below is used to perform the hypothesis testing.
L_ck.csr <- envelope(childcare_ck_ppp, Lest, # generate L-function CSR envelopes for CK
nsim = 99, # 99 simulations (as per screenshot)
rank = 1, # rank-1 global envelope
glocal = TRUE) # global+local envelope optionGenerating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
99.
Done.
plot(L_ck.csr, . - r ~ r, # plot (L_hat - r) vs r using Prof’s plotting formula
xlab = "d", ylab = "L(d)-r") # axis labels exactly as in the slide
5.9.2 Tampines planning area
5.9.2.1 Computing L-function estimation
L_tm <- Lest(childcare_tm_ppp, correction = "Ripley") # estimate L(r) for Tampines
plot(L_tm, . - r ~ r, # plot (L_hat - r) vs r as in Prof’s figure
ylab = "L(d)-r", xlab = "d(m)", # axis labels to match the slide
xlim = c(0, 1000)) # distance window 0–1000 m as shown
5.9.2.2 Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
\(H_0\) = The distribution of childcare services at Tampines are randomly distributed.
\(H_1\) = The distribution of childcare services at Tampines are not randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.
The code chunk below will be used to perform the hypothesis testing.
L_tm.csr <- envelope(childcare_tm_ppp, Lest, # CSR envelopes for L(r) in Tampines
nsim = 99, # 99 simulations
rank = 1, # rank-1 global envelope
glocal = TRUE) # global+local envelope handlingGenerating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
99.
Done.
plot(L_tm.csr, . - r ~ r, # plot (L_hat - r) vs r per Prof’s code
xlab = "d", ylab = "L(d)-r", # keep labels consistent with the slide
xlim = c(0, 500)) # show 0–500 m window exactly as given